home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 1.iso / toolbox / src / demos / OpenGL / stonehenge / TimeDate.c++ < prev    next >
C/C++ Source or Header  |  1996-11-11  |  9KB  |  327 lines
  1. /*
  2.  * (c) Copyright 1993, Silicon Graphics, Inc.
  3.  * ALL RIGHTS RESERVED
  4.  *
  5.  * Permission to use, copy, modify, and distribute this software for
  6.  * any purpose and without fee is hereby granted, provided that the above
  7.  * copyright notice appear in all copies and that both the copyright notice
  8.  * and this permission notice appear in supporting documentation, and that
  9.  * the name of Silicon Graphics, Inc. not be used in advertising
  10.  * or publicity pertaining to distribution of the software without specific,
  11.  * written prior permission.
  12.  *
  13.  * THE MATERIAL EMBODIED ON THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"
  14.  * AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,
  15.  * INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR
  16.  * FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL SILICON
  17.  * GRAPHICS, INC.  BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,
  18.  * SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY
  19.  * KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,
  20.  * LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF
  21.  * THIRD PARTIES, WHETHER OR NOT SILICON GRAPHICS, INC.  HAS BEEN
  22.  * ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON
  23.  * ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE
  24.  * POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.
  25.  *
  26.  * U.S. GOVERNMENT RESTRICTED RIGHTS LEGEND
  27.  * Use, duplication, or disclosure by the Government is subject to
  28.  * restrictions set forth in FAR 52.227.19(c)(2) or subparagraph
  29.  * (c)(1)(ii) of the Rights in Technical Data and Computer Software
  30.  * clause at DFARS 252.227-7013 and/or in similar or successor
  31.  * clauses in the FAR or the DOD or NASA FAR Supplement.
  32.  * Unpublished-- rights reserved under the copyright laws of the
  33.  * United States.  Contractor/manufacturer is Silicon Graphics,
  34.  * Inc., 2011 N.  Shoreline Blvd., Mountain View, CA 94039-7311.
  35.  *
  36.  * OpenGL(TM) is a trademark of Silicon Graphics, Inc.
  37.  */
  38. #include <GL/glu.h>
  39. #include <GL/glx.h>
  40.  
  41. #include "TimeDate.h"
  42.  
  43. #include <stdio.h>
  44. #include <stdlib.h>
  45. #include <sys/time.h>
  46.  
  47. inline float sign(float a) {return a > 0.0 ? 1. : (a < 0.0 ? -1. : 0.);}
  48.  
  49. inline float degrees(float a) {return a * 180.0 / M_PI;}
  50. inline float radians(float a) {return a * M_PI / 180.0;}
  51.  
  52. inline float arcsin(float a) {return atan2(a, sqrt(1.-a*a));}
  53.  
  54. /* Thirty days hath September... */
  55. const int mday[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  56.  
  57. TimeDate::TimeDate(int hour, int minute, int second, long usecond)
  58. {
  59.   read_time();
  60.   t.tm_hour = hour;
  61.   t.tm_min = minute;
  62.   t.tm_sec = second;
  63.   usec = usecond;
  64. }
  65.  
  66. TimeDate TimeDate::operator=(TimeDate a)
  67. {
  68.   usec = a.usec;
  69.   t = a.t;
  70.   return *this;
  71. }
  72.  
  73. TimeDate TimeDate::operator+(TimeDate a)
  74. {
  75.   TimeDate td;
  76.   td = *this + a.t;
  77.   td.usec = usec + a.usec;
  78.   return td;
  79. }
  80.  
  81. TimeDate TimeDate::operator+(struct tm t2)
  82. {
  83.   TimeDate td;
  84.   td.usec = 0;
  85.   td.t.tm_sec = t.tm_sec + t2.tm_sec;
  86.   td.t.tm_min = t.tm_min + t2.tm_min;
  87.   td.t.tm_hour = t.tm_hour + t2.tm_hour;
  88.   td.t.tm_mday = t.tm_mday + t2.tm_mday;
  89.   td.t.tm_mon = t.tm_mon + t2.tm_mon;
  90.   td.t.tm_year = t.tm_year + t2.tm_year;
  91.   td.t.tm_wday = t.tm_wday + t2.tm_wday;
  92.   td.t.tm_yday = t.tm_yday + t2.tm_yday;
  93.   return td;
  94. }
  95.  
  96. TimeDate TimeDate::operator+=(TimeDate a)
  97. {
  98.   /* This is slower than it needs to be */
  99.   *this = *this + a.t;
  100.   return *this;
  101. }
  102.  
  103. TimeDate TimeDate::operator-(TimeDate a)
  104. {
  105.   TimeDate td;
  106.   td.usec = usec - a.usec;
  107.   td.t.tm_sec = t.tm_sec - a.t.tm_sec;
  108.   td.t.tm_min = t.tm_min - a.t.tm_min;
  109.   td.t.tm_hour = t.tm_hour - a.t.tm_hour;
  110.   td.t.tm_mday = t.tm_mday - a.t.tm_mday;
  111.   td.t.tm_mon = t.tm_mon - a.t.tm_mon;
  112.   td.t.tm_year = t.tm_year - a.t.tm_year;
  113.   td.t.tm_wday = t.tm_wday - a.t.tm_wday;
  114.   td.t.tm_yday = t.tm_yday - a.t.tm_yday;
  115.   return td;
  116. }
  117.  
  118. TimeDate TimeDate::operator*(float f)
  119. {
  120.   TimeDate td;
  121.   td.usec = (long)((float)usec * f);
  122.   td.t.tm_sec = (int)((float)t.tm_sec * f);
  123.   td.t.tm_min = (int)((float)t.tm_min * f);
  124.   td.t.tm_hour = (int)((float)t.tm_hour * f);
  125.   td.t.tm_mday = (int)((float)t.tm_mday * f);
  126.   td.t.tm_mon = (int)((float)t.tm_mon * f);
  127.   td.t.tm_year = (int)((float)t.tm_year * f);
  128.   td.t.tm_wday = (int)((float)t.tm_wday * f);
  129.   td.t.tm_yday = (int)((float)t.tm_yday * f);
  130.   return td;
  131. }
  132.  
  133. TimeDate TimeDate::read_time()
  134. {
  135.   struct timeval tv;
  136.   time_t cal_time;
  137.   struct tm *tmp;
  138.  
  139.   cal_time = time(NULL);
  140.   gettimeofday(&tv);
  141.   tmp = localtime(&cal_time);
  142.  
  143.   usec = tv.tv_usec;
  144.   t = *tmp;
  145.  
  146.   return *this;
  147. }
  148.  
  149. void TimeDate::print()
  150. {
  151.   puts(asctime(&t));
  152. }
  153.  
  154.  
  155. /* The sun_direction routine comes from an awk program by Stuart Levy
  156.  * (slevy@geom.umn.edu) */
  157.  
  158. /* Day number of March 21st, roughly the vernal equinox */
  159. const float equinox = mday[0] + mday[1] + 21;
  160.  
  161. /*
  162.  *        Date     E.T. [degrees]
  163.  *        Jan 1    0.82
  164.  *        Jan 31   3.35
  165.  *        Mar 1    3.08
  166.  *        Mar 31   1.02
  167.  *        Apr 30  -0.71
  168.  *        May 30  -0.62
  169.  *        Jun 29   0.87
  170.  *        Jul 29   1.60
  171.  *        Aug 28   0.28
  172.  *        Sep 27  -2.28
  173.  *        Oct 27  -4.03
  174.  *        Nov 26  -3.15
  175.  *        Dec 26   0.19
  176.  */
  177. const float et[] = {
  178.   .82, 3.35, 3.08, 1.02, -.71, -.62, .87,
  179.   1.60, .28, -2.28, -4.03, -3.15, .19, 2.02};
  180.  
  181. const float DEG = 180. / M_PI;
  182.  
  183. Point TimeDate::sun_direction(float lon, float lat)
  184. {
  185.   float yearday, eti, etoffset, ET, HA, LON, decl, altitude, azimuth;
  186.   int etindex;
  187.   Point dir;
  188.   long tmp_usec;
  189.   int m;
  190.  
  191.   tmp_usec = usec;
  192.   usec = 0;
  193.   *this = correct_smaller();
  194.  
  195.   /*
  196.    *   hour angle of the sun (time since local noon):
  197.    *       HA = (time - noon) * 15 degrees/hour
  198.    *          + (observer''s east longitude - long. of time zone''s 
  199.    *                   central meridian
  200.    *                           POSITIVE east, NEGATIVE west)
  201.    *          - (15 degrees if time above was daylight savings time, else 0)
  202.    *          - "equation of time" (depends on date, correcting the 
  203.    *           hour angle for
  204.    *               the earth''s elliptical orbit and the slope of the ecliptic)
  205.    */
  206.   yearday = (float)t.tm_mday + 
  207.     ((float)t.tm_hour + (float)t.tm_min / 60.0) / 24.0;
  208.   for (m = 0; m < t.tm_mon; m++) yearday += mday[m];
  209.  
  210.   /*
  211.    * Index into equation-of-time table, tabulated at 30-day intervals
  212.    * We''ve added an extra entry at the end of the et table, corresponding
  213.    * to January 25th of the following year, to make interpolation work
  214.    * throughout the year.
  215.    * Use simple linear interpolation.
  216.    */
  217.   eti = (yearday - 1.) / 30.;
  218.   etoffset = eti - trunc(eti);
  219.   etindex = (int)trunc(eti) + 1;
  220.   ET = et[etindex - 1] + etoffset*(et[etindex+1 - 1] - et[etindex - 1]);
  221.   /* The 90. puts us in the Central time zone */
  222.   HA = ((float)t.tm_hour + (float)t.tm_min/60. - 12.)*15. +
  223.     lon + 90. - ET;
  224.  
  225.   /*
  226.    *    Sun''s declination:
  227.    *        sun''s longitude = (date - March 21) * 360 degrees / 365.25 days
  228.    *                [This ignores the earth''s elliptical orbit...]
  229.    */
  230.   LON = (yearday - equinox) * 360. / 365.25;
  231.   decl = DEG * arcsin( sin(LON/DEG) * sin(23.4/DEG) );
  232.  
  233.   /* 
  234.    *    Then a spherical triangle calculation to convert the Sun''s
  235.    *    hour angle and declination, and the observer''s latitude,
  236.    *    to give the Sun''s altitude and azimuth (angle east from north).
  237.    */
  238.   altitude = DEG * arcsin( sin(decl/DEG)*sin(lat/DEG) 
  239.               + cos(decl/DEG)*cos(lat/DEG)*cos(HA/DEG) );
  240.   azimuth = DEG * atan2( -cos(decl/DEG)*sin(HA/DEG), 
  241.             sin(decl/DEG)*cos(lat/DEG) 
  242.             - cos(decl/DEG)*cos(HA/DEG)*sin(lat/DEG) );
  243.  
  244.   /*
  245.   printf("On %d/%d at %d:%02d, lat %g, lon %g\n", 
  246.      t.tm_mon + 1, t.tm_mday + 1, t.tm_hour, t.tm_min, lat, lon);
  247.   printf("HA %.1f ET %.1f Sun lon %.1f decl %.1f alt %.1f az %.1f\n", 
  248.      HA,ET,LON,decl,altitude,azimuth);
  249.   */
  250.   
  251.   dir.pt[2] = sin(radians(altitude));
  252.   dir.pt[0] = cos(radians(azimuth));
  253.   dir.pt[1] = sin(radians(azimuth));
  254.   dir.pt[3] = 0;
  255.   dir.unitize();
  256.   // dir.print();
  257.  
  258.   usec = tmp_usec;
  259.  
  260.   return dir;
  261. }
  262.  
  263. TimeDate TimeDate::correct_bigger()
  264. {
  265.   TimeDate td;
  266.   int days;
  267.  
  268.   td.t = t;
  269.   td.usec = usec;
  270.  
  271.   td.t.tm_sec += (int)(td.usec / 1000000);
  272.   td.usec %= 1000000;
  273.   td.t.tm_min += td.t.tm_sec / 60;
  274.   td.t.tm_sec %= 60;
  275.   td.t.tm_hour += td.t.tm_min / 60;
  276.   td.t.tm_min %= 60;
  277.  
  278.   if (td.t.tm_hour > 23) {
  279.     days = td.t.tm_hour / 24;
  280.     td.t.tm_hour %= 24;
  281.     td.t.tm_mday += days;
  282.     td.t.tm_wday = (td.t.tm_wday + days) % 7;
  283.     td.t.tm_yday = (td.t.tm_yday + days) % 365;
  284.     while (td.t.tm_mday > mday[td.t.tm_mon]) {
  285.       td.t.tm_mday -= mday[td.t.tm_mon++];
  286.       if (td.t.tm_mon >= 12) {
  287.     td.t.tm_year += td.t.tm_mon / 12;
  288.     td.t.tm_mon %= 12;
  289.       }
  290.     }
  291.   }
  292.   return td;
  293. }
  294.  
  295. TimeDate TimeDate::correct_smaller()
  296. {
  297.   TimeDate td;
  298.  
  299.   td.t = t;
  300.   td.usec = usec;
  301.  
  302.   while (td.usec < 0) {
  303.     td.t.tm_sec--;
  304.     td.usec += 1000000;
  305.   }
  306.   while (td.t.tm_sec < 0) {
  307.     td.t.tm_min--;
  308.     td.t.tm_sec += 60;
  309.   }
  310.   while (td.t.tm_min < 0) {
  311.     td.t.tm_hour--;
  312.     td.t.tm_min += 60;
  313.   }
  314.   while (td.t.tm_hour < 0) {
  315.     td.t.tm_mday--;
  316.     if (td.t.tm_wday) td.t.tm_wday--;
  317.     else td.t.tm_wday = 6;
  318.     td.t.tm_yday--;
  319.     td.t.tm_hour += 24;
  320.   }
  321.   if (td.t.tm_mon < 0 || td.t.tm_year < 0 || td.t.tm_yday < 0)  {
  322.     fprintf(stderr, "Warning:  day < 0 in TimeDate.c++.\n");
  323.     td.t.tm_mon = td.t.tm_year = td.t.tm_yday = 0;
  324.   }
  325.   return td.correct_bigger();
  326. }
  327.